Parallel wave-based analog computing using metagratings

Abstract Wave-based signal processing has witnessed a significant expansion of interest in a variety of science and engineering disciplines, as it provides new opportunities for achieving high-speed and low-power operations. Although flat optics desires integrable components to perform multiple missions, yet, the current wave-based computational metasurfaces can engineer only the spatial content of the input signal where the processed signal obeys the traditional version of Snell’s law. In this paper, we propose a multi-functional metagrating to modulate both spatial and angular properties of the input signal whereby both symmetric and asymmetric optical transfer functions are realized using high-order space harmonics. The performance of the designed compound metallic grating is validated through several investigations where closed-form expressions are suggested to extract the phase and amplitude information of the diffractive modes. Several illustrative examples are demonstrated to show that the proposed metagrating allows for simultaneous parallel analog computing tasks such as first- and second-order spatial differentiation through a single multichannel structured surface. It is anticipated that the designed platform brings a new twist to the field of optical signal processing and opens up large perspectives for simple integrated image processing systems.


Introduction
For a few decades, digital processors have been widely used to execute computational tasks, as an alternative to analog mechanical and electrical computers. Despite their reliability and high-speed operation, digital processors suffer from high-power consumption, expensive analogto-digital conversion, and sharp performance degradation at high frequencies, leading to large limitations even for performing simple computing tasks such as differentiation or integration, equation solving, matrix inversion, edge detection, and image processing [1,2]. With the advent of metamaterials and metasurfaces, spatial analog optical computing resurfaced, finding important applications as compact solutions for high speed, high throughput image processing and parallel computing. Since the seminal proposal of Silva et al., [3] wave-based analog computing has witnessed rapid progress, with the demonstration of optical spatial differentiators [4][5][6][7][8][9][10][11][12][13], integrators [2,8,12,14], equation solvers [15][16][17], spatiotemporal computing [9,18,19], and wave-based neuromorphic computing [20][21][22][23]. Among them, the Green's function (GF) method in which a specific-purpose computing operation is directly realized in the real space, without transforming back and forth from the spatial to the spectral domain [14], affords compactness and avoids possible challenges in error propagation and alignment issues. The applicability of the GF method to execute signal processing has been verified in a series of proposals via spin hall effect of light [24], disordered and complex scattering system [7,25,26], layered structures [18,27], topological insulators [28,29], plasmonic arrays [5], bianisotropic metasurfaces [10,17], and so on. Nevertheless, prior GF-based studies still face two different challenges: (i) parallel realization of mathematical operators has been only addressed by using bulky structures [7,16] and array of subwavelength meta-atoms with complex geometries [30,31] and thus, they are still subject to implementation difficulties arising from high fabrication precision demands; (ii) although reflective optical processing for normal incidences is a good alternative for complex oblique illumination setups, it still needs additional optical components to separate the processed signal from the input one [10]. Further efforts to tackle these barriers must be accompanied with the use of more powerful architectures to implement spatial optical signal processing.
Passive metasurfaces can ensure highly-efficient wavefront molding by leveraging nonlocal effects stemming from the excitation of evanescent [32,33] or leaky modes [34], but this typically leads to complex design requirements and a need for deeply subwavelength fabrication resolution [35][36][37][38][39][40][41][42][43][44]. As an alternative solution, recent efforts have shown that metagratings composed of non-subwavelength periodic patterns enable highly complex diffraction scenarios of significant practical interest, with high efficiency [45,46]. The core idea of metagratings is based on Floquet-Bloch (FB) theory, remarking that when a plane wave impinges on a periodic structure, a discrete set of diffracted waves can be generated, some of them propagating and others being evanescent [47,48]. The number of propagating and evanescent waves is determined by the period of the structure and the angle of incidence. The structure is engineered at the scale of the wavelength so as to suppress the unwanted space harmonics and reroute the incident power towards a desired nonspecular channel. Therefore, such design potentially relaxes some of the fabrication challenge of metasurfaces, as it does not require precise lithography techniques [45,49]. Metagratings have also the ability to provide multiple arbitrarily-oriented space channels for creating multimission surfaces [50]. Although several reports have examined metagratings from different points of view, the potential application of these structures for performing optical analog computation, which requires studying the phase information of the space harmonics, has not been unveiled, yet.
In this study, we reveal that a suitably engineered all-metallic metagrating can open multiple non-specular channels for parallel implementation of symmetric and asymmetric optical transfer functions at both normal and oblique illuminations. Closed-form expressions are presented to predict the reflection phase and amplitude information of the space channels opened by the designed compound metallic grating. The application of the designed metagrating for accomplishing first-and second-order differentiation as well as detecting sharp edges of an input image is investigated through several illustrative demonstrations.

Fundamental theory of compound metallic grating
The main idea of this paper is graphically depicted in Figure 1a, where a single metagrating is responsible for realizing parallel mathematical operations on the input signals coming from different directions. The symmetric geometry originated from the presence of a single groove in each period does not allow realizing the asymmetry required for performing odd-order transfer functions. Thus, we consider a two-dimensional (2D) sparse array ( ∕ y = 0) composed of two grooves in each period of a metallic medium filling the half-space z < 0, as shown in Figure 1b  , respectively. The whole structure is surrounded by a medium with dielectric constant r0 ( . The position and dimensions of the grooves can be set at will and the periodicity along x axis is indicated by L x . The metagrating is illuminated by a general oblique 2D beam profile with the angle in and the transverse-magnetic (TM) polarization. By leveraging the superposition principle and spectral decomposition, the input f in and output f out optical signals can be expanded based on an infinite set of plane waves, all traveling in different directions with wavenumbers k x : Here, W denotes the beamwidth of the input field and the harmonic time dependency e j t is omitted. The spatial frequency content of the incident field is represented by the various plane wave amplitudes around k x = k 0 sin( ), that together form the input signal. The compound metallic grating, depicted in Figure 1c, interacts differently with each of these plane wave components. A Rayleigh expansion can be performed to write the total fields at the upper half-space: in which, the subscript m corresponds to the order of space harmonics, and Here, k x0 = k 0 n 0 sin( i ) and Y m = k zm ∕ 0 n 2 0 denote respectively the x-directed wavenumber the admittance of the mth TM-polarized diffractive mode in the upper half-space. After applying the proper boundary conditions and solving the related equations, the specular and nonspecular reflection coefficients can be expressed as: in which, A, A m , and B are complex numbers. A detailed derivation of the mathematical expressions of these coefficients can be found in Supplementary Information A. It should be noted that Eqs. (4a) and (4b) include both phase and amplitude information of the diffractive modes when the metagrating is excited by an oblique plane wave with an arbitrary angle of incidence.
manipulations on Eq. (5), the unlocking condition for each diffractive mode turns into:  6), disclosing the best choice for grating periodicity corresponding to the desired angle of incidence. As designed, only three FB modes m = 0, ±1 will be propagating for the input beam profile, and the rest will be evanescent, as long as the periodicity and the incident wave angle satisfy Hereafter, we intend to show how the multiple channels provided by the metagrating can be exploited for performing analog signal processing and be used for realizing different functionalities at the same time.

Single-operator metasurface
We adjust the periodicity of the metagrating so as to provide three active channels (see Figure 2a). Accordingly, we set the periodicity of the designed metagrating as L x = 1.3 0 (f 0 = 1 THz), which orients the channels along = ±40 • directions at f = 1.2 THz. Once the input and output channels are determined, the width and height of the contributing grooves are optimized so that the angular dispersion of the scattering coefficient between these channels in Eqs. (4a) and (4b), emulates the k xdependency of desired transfer functions. For instance, we can achieve 1st-and 2nd-order spatial differentiation of the input signal exciting the metasurface from the ith channel and exiting at the jth channel, provided that S ji ( ) = jk 0 (sin( ) − sin( inc )) and S ji ( ) = −k 2 0 (sin( ) − sin( inc )) 2 , respectively. Here, our goal is to design a periodic surface that applies 1st-order spatial differentiation operation on the input signal traveling from port 2 to port 1 (see Figure 2b). A comprehensive parametric study based on the theoretical representation of Eqs. (4a) and (4b), has been carried out to find the optimum parameters of the metagrating and minimize the following error function: Here,H( ) indicates the transfer function of choice, i and j refer to the output and input ports, respectively, and the above summation is calculated over q discrete angles in the vicinity of the incident wave angle. The angular spectra of |S 12 ( )| and ∠S 12 ( ) are plotted in To evaluate the performance of our single-operator metasurface differentiator, the Gaussian signal of Figure 2e is utilized to launch the metagrating from port 2 and the calculated output field is shown in the same figure. The result is compared with the exact response indicating that the output field at port 1 is indeed the response expected from a first-order derivative operation, a 2% error, as defined by Eq. (7). The error value of each computing task is given in the caption of the related figure. Indeed, the designed metagrating successfully deflects the first-order derivative of the normally incident signals into = −40 • direction without using any additional beam splitting sub-block. Conversely, previous metasurfaces based on the GF method cannot provide a reciprocal solution for asymmetric on-axis transfer functions at the reflection side [10]. Our structure therefore relaxes vexing complexities of previous designs that needed oblique illumination setups and beam splitting devices [10].

Multi-operator metasurface
Up to now, the proposed metagratings realize only a single processing operation on a single channel. Hereafter, we intend to design metagratings which successfully create multiple channels, each of which enables different processing or scattering functionalities. To this aim, the optimization procedure of Eq. (7) is simultaneously accomplished for different scattering parameters, e.g., S ij and S uv . To exemplify how the operation frequency can be adjusted, we decided to change the periodicity of the structure to L x = 1.22 0 , so that the orientation of m = ±1 channels points in the = 40 • direction for normal excitations at a higher frequency, namely f = 1.3 THz.
Both specular and nonspecular channels can be involved and we carry out the numerical optimizations for different computing scenarios displayed in Figures 3a, 4a, and 5. In the first demonstration, the metagrating is parametrically adjusted so that the S 13 ( ) (specular) and S 22 ( ) (specular) coefficients implement first-and second-order differentiation operations, respectively (see Figure 3a). The synthesized transfer functions are plotted in Figure 3b and c. A parabolic shape is acquired for |S 22  Eqs. (1a), (1b) and the achieved scattering coefficients, the output fields leaving the metagrating from ports 2 and 3 are calculated and presented in Figure 3d and e, respectively. As seen, the output signals at ports 2 and 3 are nothing but the first-and second-order derivative of the input fields at ports 1 and 2, respectively. The obtained results indicate that a suitably-designed metagrating successfully provides multi-input multi-output channels, each of which, enables different processing functionalities. Our purpose in the second example is to design a single-input multi-output metagrating manifesting the angular dispersion of jk 0 sin (first-order spatial differentiation) in its S 21 ( ) (nonspecular) and S 31 ( ) (specular) scattering coefficients (see Figure 4a). The transfer functions realized by the optimized metagrating are shown in Figure 4b and c from which, one can immediately deduce that both S 21 ( ) and S 31 ( ) parameters expose a linear trend and an asymmetric 180 • phase jump near their corresponding angles, consistent with the transfer function of the 1st-order differentiation operation. The structure has one input field but the output fields at ports 2 and 3 are examined by input signals with profiles of sin( x) and sinc( x), respectively. The input and output fields displayed in Figure 4d and e verify that this single-input multi-output processing mission is perfectly accomplished by the designed metagrating. For a quantitative comparison, the exact output responses are also illustrated. The demonstrations presented in Figures 3  and 4 show different aspects of parallel processing ability of the designed metagratings. Another example of parallel computing can be found in Supplementary Information D in which the first-order spatial differentiation can be applied on two distinct input signals at the same time. As a final demonstration example, we form a multi-operator  metagrating performing second-order spatial differentiation for the input signals coming from either port 1 or 2 (see Figure 5a). Figure 5b and c demonstrate the angledependent amplitude and phase of the transfer function of the designed metagrating. All indicator features of the second-order differentiation operation are observable in S 22 ( ) (specular) and S 11 ( ) (nonspecular) around = 0 • and = 40 • , respectively. Figure 5d and e depicts the input fields, exp(− x 2 ) and exp(− x 4 ), and the corresponding output signals, confirming that the designed multi-operator metagrating is capable of independently processing two input signals.
To further validate the performance of the proposed metagratings, full-wave simulations have been carried out through COMSOL Multiphysics. The configuration is shown in Figure 6g in which two TM-polarized Gaussian-shape beams illuminate the metagrating of Figure 6a from two different directions = 0 • and = 40 • (see Figure 6a and d). The medium surrounding the designed metagratings is filled by air and the boundary conditions are selected as perfect match layer (PML). Figure 6b and d demonstrate the scattered fields for each of normal and oblique illuminations. A cut-line of the scattered fields in each case is plotted in Figure 6c and f indicating that the output signals successfully obey the exact version of the second-order spatial derivative of the Gaussian-shape beams. Thus, the full-wave simulations verify well our idea and analytical model, confirming the relevance of high-order Floquet modes to perform analog signal processing. To investigate the possible errors originating from the fabrication tolerances, we made some geometrical variations on the structure and the corresponding results are presented in Supplementary Information D. Finally, a concrete 1D edge detection application is evaluated by using the metagrating differentiator of Figure 2b, where the logo of "EPFL" university is utilized as the input image (see Figure 7a). Due to the 1D nature of the proposed system, the reflected image displayed in Figure 7b highlights all edges of the incident image along the vertical direction. Supplementary Information E depicts another interesting application in the field of image processing in which the metagrating applies a denoising filter on the input image. It should be noted that the idea can be extended to 2D scenarios, including differentiation and edge detection, the reflection coefficients of two-dimensional metagratings should be extracted,  allowing us to manipulate the information of both spatial frequencies kx and ky [53][54][55]. Besides, the proposed multi-operator metagrating can also be an important key in parallel image processing applications, including THz imaging. With the ability to penetrate many materials, terahertz waves have intriguing potential in revealing hidden features [56]. In order to extract the important geometric features of objects, the edge detection method can be directly employed through different THz imaging systems [57]. As a result, THz-based edge detectors like those we proposed in the submitted paper, may be a key component in THz imaging/image processing systems. As a promising outlook, by tuning the interference between the multipoles of all-dielectric metagratings, one can also move this kind of signal processing platforms to optical frequencies.

Conclusions
In summary, we exploited high-order spatial harmonics of a multi-functional angle-multiplexed metagrating to realize various scattering and signal processing functionalities at the same time. Analytical expressions were presented for engineering the phase and amplitude information of each spatial harmonic. We demonstrated metagrating configurations that implement diverse optical computing functions such as first-and second-order differentiation operations, and several scattering functionalities such as retro-reflection when excited from different channels. Full-wave finite-element simulations confirmed our theory. Using the proposed angle-multiplexed metagratings, we completely escape the problem of realizing asymmetric optical transfer functions for normal illuminations, avoiding bulky additional splitting blocks to separate the routes of the incident and reflected signals, and allowing for much more compact configurations compared with previous attempts. Besides, the sparsity of our wavelength-scale array significantly relaxes fabrication tolerance, which is a crucial advantage to move the concept towards optical frequencies. This new class of multifunctional processors may find great potential applications in integrated photonic devices and imaging systems that process optical signals coming from different directions at the speed of light, as they reflect onto structured surfaces.
Author contributions: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.